      PROGRAM SOLPOW(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
C
C                   *** SOLAR POWER ARRAY MODEL ***
C
      COMMON /FTAC/C(25,25),CAP(25),Q(25),T(25)
      COMMON /ORBC/A,E,RA,PHI,OMEGA
      DIMENSION XA(3,10),AA(10),VNA(3,10),XE(3,520),AE(520),VNE(3,520),
     *          XS(3,520),AS(520),VNS(3,520),XEO(3),XO(3)
      DIMENSION VNB(3,10)
      DIMENSION TP(100),TH(100),TPE(100),THE(100)
      DIMENSION FES(520),XD(3),AD(1),VD(3)
      DIMENSION QD(500),QR(500),QT(500),TA(500),ET(500),PS(500),TM(500),
     *          TITLE(24)
      DIMENSION XP(520),YP(520),ZP(520),XR(500),YR(500),ZR(500),
     *          XSOL(2),YSOL(2),ZSOL(2)
      DIMENSION PLOAD(50),TML(50),QARRAY(500)
C
      NAMELIST /ARRAY/T,AREA,E1,E2,A1,A2,WSIL,WPOL,WADH,ALIFE
C                    T=TEMPERATURE OF ARRAY (R)
C                    AREA =AREA OF ARRAY (FT2)
C                    E1=EMITTANCE OF CELL SIDE
C                    E2=EMITTANCE OF BACK SIDE
C                    A1=SOLAR ABSORPTANCE OF CELL SIDE
C                    A2=SOLAR ABSORPTANCE OF BACK SIDE
C                    WSIL=WEIGHT OF FUSED SILICA (LB)
C                    WPOL=WEIGHT OF POLYIMIDE SUBSTRATE (LB)
C                    WADH=WEIGHT OF ADHESIVE (LB)
C                    ALIFE=FRACTION OF ARRAY DESIGN LIFE
C
      NAMELIST /LOAD/PLOAD,TML,NL
C                    PLOAD=LOAD POWER LEVEL (KW)
C                    TML  =TIME (HR)
C                    NL   =NUMBER OF POINTS
C
      NAMELIST /INPUT/DATE,A,E,RA,PHI,OMEGA,NA,NE,NS,NT,KRERN
C                    DATE=DATE OF THE YEAR (U.T.)
C                               EXAMPLE, 6.215=12 NOON ON JUNE 21 (U.T.)
C                    A   =SEMIMAJOR AXIS (FT)
C                    E   =ECCENTRICITY
C                    RA  =RIGHT ASCENSION
C                    PHI=INCLINATION
C                    OMEGA=ARGUMENT OF PERIGEE
C                    NA  =NUMBER OF ARRAY ELEMENTS
C                    NE  =NUMBER OF EARTH LATITUDE DIVISIONS
C                    NS  =NUMBER OF SOLAR LATITUDE DIVISIONS
C                    NT  =NUMBER OF TIME STEPS PER ORBIT
C                    KRERN=0 FOR STOP, 1 FOR NEW INPUT, 2 FOR ALL NEW CASE
C
C
C
      DATA C,CAP,Q,JRKS/675*0.0,1/
      DATA SIGMA/1.714E-9/
C
      DATA ZSOL/2*0.0/
C
C
C                  EARTH ALBEDO,EMITTANCE
       DATA BEARTH,EEARTH/0.30,0.86/
C
C                  TEMPERATURE OF EARTH,SUN,COSMIC BACKGROUND RADIATION
       DATA T/0.,   479.1, 10467.0, 4.97,   21*0.0/
C           EARTH ORBIT PARAMETERS
      DATA XEO,XO/6*0.0/
      DATA AEARTH,EE,WS,WE/4.90745E11,0.016718,4.387E30,1.317E25/
C           EARTH,SOLAR RADII
      DATA RE,RS/2.08982E7,2.28228E9/
C
    2    CONTINUE
      READ(5,ARRAY)
    1    CONTINUE
      READ(5,LOAD)
      READ(5,INPUT)
      READ(5,3) TITLE
    3 FORMAT(8A10)
C
C                      DETERMINE ORBITAL EPHEMERIDES
C                       ...EARTH
                            CALL EPHEM(AEARTH,EE,WS,PE,TPE,THE)
C                       ...SOLAR ARRAY
                            CALL EPHEM(A,E,WE,P,TP,TH)
C
C                      LOCATION OF EARTH IN HELIOCENTRIC COORDINATES
C                      ... CENTER
                            CALL EARTH(DATE,XEO(1),XEO(2),TPE,THE,PE)
C                      ... SURFACE ELEMENTS SEEING ORBITAL TRACK
                                 DPHI=ACOS(RE/(A*(1.+E)))
                                 PHIMIN=PI(0.5)-DPHI
                                 PHIMAX=PI(0.5)+DPHI
                         CALL SPHERE2(XEO,PHI,RA,RE,XE,AE,VNE,
     *                                PHIMIN,PHIMAX,NE,NEA)
C                      SOLAR BETA ANGLE
                         SBETA=BETA(XEO,PHI,RA)
C                      SOLAR MODEL
                         CALL SPHERE(XO,RS,XS,AS,VNS,0.,PI(1.),NS,NSA)
C                      ELEMENT SHAPE FACTORS FOR ALBEDO CALCULATION
      SIGAE=0.0
         DO 10 J=1,NEA
                XD(1)=XE(1,J)   $      VD(1)=VNE(1,J)
                XD(2)=XE(2,J)   $      VD(2)=VNE(2,J)
                XD(3)=XE(3,J)   $      VD(3)=VNE(3,J)
                              AD(1)=AE(J)
                SIGAE=SIGAE+AD(1)
                FES(J)=SHFAC(XD,AD,VD,1,XS,AS,VNS,NSA,XD,AD,VD,1)
   10    CONTINUE
C
C
          DT=P/NT
          J=1
          K=NT+1
      JLOAD=2
          TIME=0.0
   95 CONTINUE
C                      LOCATION OF ARRAY IN HELIOCENTRIC COORDINATES
                CALL ORBIT(XEO(1),XEO(2),TIME,X,Y,Z,TP,TH,P)
                            XA(1,1)=X  $  XA(2,1)=Y  $  XA(3,1)=Z
                             XR(J)=X  $  YR(J)=Y  $  ZR(J)=Z
                            AA(1)=1.0
                             R=SQRT(X*X+Y*Y+Z*Z)
                            VNA(1,1)=0.-X/R
                            VNA(2,1)=0.-Y/R
                            VNA(3,1)=0.-Z/R
                            VNB(1,1)=0.-VNA(1,1)
                            VNB(2,1)=0.-VNA(2,1)
                            VNB(3,1)=0.-VNA(3,1)
C
C                      DIRECT SOLAR RADIATION
      FSOL=SHFAC(XA,AA,VNA,NA,XS,AS,VNS,NSA,XE,AE,VNE,NEA)
             QD(J)=FSOL*SIGMA*AREA*T(3)**4.
C                      REFLECTED RADIATION
      IF(QD(J).GT.0.00)  GO TO 20
               QR(J)=0.0
               Q2=0.0
                       GO TO 30
   20 CONTINUE
C                       ... CELL SIDE
      QR(J)=FALBEDO(XA,AA,VNA,NA,XE,AE,VNE,FES,NEA)*SIGMA*AREA*BEARTH
     *                 *T(3)**4.
C                       ... BACK SIDE
      Q2   =FALBEDO(XA,AA,VNB,NA,XE,AE,VNE,FES,NEA)*SIGMA*AREA*BEARTH
     *              *A2*T(3)**4.
   30 CONTINUE
C
C                      TERRRESTIAL RADIATION
      FEA=SHFAC(XA,AA,VNA,NA,XE,AE,VNE,NEA,XD,AD,VD,1)
      FEB=SHFAC(XA,AA,VNB,NA,XE,AE,VNE,NEA,XD,AD,VD,1)
   96 CONTINUE
C
      C(1,2)=CRAD(E1,EEARTH,AREA,SIGAE,FEA)+
     *       CRAD(E2,EEARTH,AREA,SIGAE,FEB)
C
C                      COSMIC BACKGROUND RADIATION
      C(1,4)=CRAD(E1,1.,AREA,1.,1.-FEA-FSOL)+
     *      CRAD(E2,1.,AREA,1.,1.-FEB)
C
C                      LOAD POWER
      PS(J)=YRDR(TIME,TML,PLOAD,JLOAD,NL)
C                      ARRAY EFFICIENCY
                    ET(J)=ETA(T(1),ALIFE)
C                     MAXIMUM ARRAY GENERATION RATE
              QARRAY(J)=(QD(J)+QR(J))*ET(J)/3412.76
C                      ACTUAL ARRAY GENERATION RATE
                     IF(PS(J).LT.QARRAY(J))
     *        QARRAY(J)=PS(J)
C                      EFFECTIVE ARRAY HEATING RATE
              Q(1)=(QD(J)+QR(J))*A1+Q2-QARRAY(J)*3412.76
     *             -SIGMA*E1*AREA*FSOL*T(1)**4.
C                      HEAT CAPACITY OF ARRAY
      CAP(1)=WSIL*CPGLAS(1)+WPOL*CPKAP(1)+WADH*0.35
C
                    IF(JRKS.GT.1) GO TO 90
                  TM(J)=TIME
              QD(J)=QD(J)/3412.76
              QR(J)=QR(J)/3412.76
              QT(J)=QD(J)+QR(J)
                    ET(J)=0.0
                     IF(QT(J).NE.0.0)
     *              ET(J)=QARRAY(J)*100./QT(J)
                  TA(J)=T(1)-459.67
                     IF(J.EQ.K) GO TO 100
                     J=J+1
   90    CONTINUE
                     CALL STDTM(DT,TIME,JRKS)
                    IF(JRKS.EQ.1) GO TO 96
                    IF(JRKS.EQ.3) GO TO 96
                    IF(JRKS.GT.1) GO TO 95
  100                   CONTINUE
C
C                      OUTPUT
      WRITE(6,200) TITLE
  200 FORMAT(1H1,8A10/1H ,8A10/1H ,8A10)
      WRITE(6,204) SBETA
  204 FORMAT(1H0,32X,"-BETA=",F6.2," DEG-")
      WRITE(6,201) NEA,NSA
  201 FORMAT(1H0,20X,"NUMBER OF EARTH ELEMENTS=",I4/1H ,20X,"NUMBER OF S
     *OLAR ELEMENTS=",I4)
      WRITE(6,202)
  202 FORMAT(1H0,23X,"SOLAR IRRADIATION"/1H ,7X,"TIME",5X,"DIRECT  REFLE
     *CTED      TOTAL",7X,"TEMP",4X,"PWR GEN",6X,"EFFIC"/
     *1H ,8X,"HR",9X,
     *"KW",9X,"KW",9X,"KW",10X,"F",9X,"KW",8X,"PCT")
      WRITE(6,203) (TM(J),QD(J),QR(J),QT(J),TA(J),QARRAY(J),ET(J),
     *              J=1,K)
  203 FORMAT(1H ,7F11.4)
C                      EARTH SURFACE ARRAYS
      DO 50 J=1,NEA
          XP(J)=XE(1,J)
          YP(J)=XE(2,J)
          ZP(J)=XE(3,J)
   50 CONTINUE
C                      SOLAR VECTOR END POINTS
          R=SQRT(XEO(1)*XEO(1)+XEO(2)*XEO(2))
                   XSOL(1)=XEO(1)
                   XSOL(2)=1.5*RE*(0.-XSOL(1)/R)+XSOL(1)
                   YSOL(1)=XEO(2)
                   YSOL(2)=1.5*RE*(0.-YSOL(1)/R)+YSOL(1)
C
C
C
                           IF(KRERN-1) 300,1,2
  300                         CONTINUE
                                 STOP
                                  END
      REAL FUNCTION BETA(RE,PHI,RA)
C               SUBROUTINE FOR DETERMINATION OF SOLAR BETA ANGLE (DEG)
      DIMENSION RE(3)
C
C                      ORBITAL NORMAL VECTOR IN ORBIT
C                      COORDINATES
      XA=0.0
      YA=0.0
      ZA=1.0
C                      ROTATE FOR INCLINATION
      CALL RTOP(YA,ZA,R,THETA)
               THETA=THETA+PHI
      CALL PTOR(YA,ZA,R,THETA)
C                      ROTATE FOR RIGHT ASCENSION
      CALL RTOP(XA,YA,R,THETA)
               THETA=THETA+RA
      CALL PTOR(XA,YA,R,THETA)
C                      ROTATE FOR EARTH'S INCLINATION
      CALL RTOP(YA,ZA,R,THETA)
               THETA=THETA-0.409274
      CALL PTOR(YA,ZA,R,THETA)
C
       BET=ACOS((RE(1)*XA+RE(2)*YA)/SQRT(RE(1)*RE(1)+RE(2)*RE(2)))
     *     -PI(0.5)
       BETA=BET*180./PI(1.)
                       RETURN
                        END
      SUBROUTINE ORBIT(XO,YO,TIME,X,Y,Z,TP,TH,P)
C                        SUBROUTINE FOR HELIOCENTRIC
C                        COORDINATES OF AN EARTH SATELLITE
C
C               XO,YO=LOCATION OF EARTH (FT)
C               TIME=TIME FROM PERIGEE (HR)
C               X,Y,Z=SATELLITE COORDINATES (FT)
C               TP,TH=POSITION ARRAYS FROM EPHEM
C               P=   ORBITAL PERIOD (HR)
C
      COMMON /ORBC/A,E,RA,PHI,OMEGA
      DIMENSION TP(100),TH(100)
      DATA J,N/2,100/
C
C                   LOCATION IN ORBIT COORDINATES
      THETA=YRDR(TIME/P,TP,TH,J,N)
      R=A*(1.-E*E)/(1.+E*COS(THETA))
C                   ROTATE FOR ARGUMENT OF PERIGEE
      THETA =THETA+OMEGA
        CALL PTOR(X,Y,R,THETA)  $ Z=0.0
C                   ROTATE FOR INCLINATION
      CALL RTOP(Y,Z,R,THETA)
          THETA=THETA+PHI
        CALL PTOR(Y,Z,R,THETA)
C                   ROTATE FOR RIGHT ASCENSION
        CALL RTOP(X,Y,R,THETA)
          THETA=THETA+RA
        CALL PTOR(X,Y,R,THETA)
C                   ROTATE FOR EARTH'S INCLINATION
        CALL RTOP(Y,Z,R,THETA)
           THETA=THETA-0.409274
        CALL PTOR(Y,Z,R,THETA)
C                   TRANSLATE FOR EARTH'S POSITION
        X=X+XO
        Y=Y+YO
                    RETURN
                     END
      SUBROUTINE EARTH(DATE,X,Y,TP,TH,P)
C               SUBROUTINE FOR HELIOCENTRIC POSITION
C               OF EARTH ON A GIVEN DATE
C
C              DATE= DATE OF YEAR (U.T.)
C              X,Y=HELIOCENTRIC COORDINATES (FT)
C              TP,TH=POSITION ARRAYS FORM EPHEM
C              P   =ORBITAL PERIOD (HRS)
C
      DIMENSION TP(100),TH(100)
      DATA J,N/2,100/
      DATA A,E/4.90744E11,0.016718/
C
      JD=DATE
               GO TO (1,2,3,4,5,6,7,8,9,10,11,12), JD
    1 XN=100.*(DATE-1)-1.           $ GO TO 20
    2 XN=100.*(DATE-2)+30.          $ GO TO 20
    3 XN=100.*(DATE-3)+58.242       $ GO TO 20
    4 XN=100.*(DATE-4)+89.242       $ GO TO 20
    5 XN=100.*(DATE-5)+119.242      $ GO TO 20
    6 XN=100.*(DATE-6)+150.242      $ GO TO 20
    7 XN=100.*(DATE-7)+180.242      $ GO TO 20
    8 XN=100.*(DATE-8)+211.242      $ GO TO 20
    9 XN=100.*(DATE-9)+242.242      $ GO TO 20
   10 XN=100.*(DATE-10)+272.242     $ GO TO 20
   11 XN=100.*(DATE-11)+303.242     $ GO TO 20
   12 XN=100.*(DATE-12)+333.242     $ GO TO 20
C
   20 CONTINUE
       XN=24.*XN
        DTNODE=YRDR(4.49854199,TH,TP,J,N)*P-6381.091333
        XN=XN+DTNODE
              IF(XN.GT.P) XN=XN-P
              IF(XN.LT.0.0) XN=XN+P
      THETA=YRDR(XN/P,TP,TH,J,N)+1.78464317
        R=A*(1.-E*E)/(1.+E*COS(THETA))
      CALL PTOR(X,Y,R,THETA)
                            RETURN
                             END
      SUBROUTINE EPHEM(A,E,W,P,TP,THETA)
C                   SUBROUTINE FOR ORBITAL POSITION AS
C                   A FUNCTION OF TIME
C
C            A=SEMIMAJOR AXIS (FT)
C            E=ECCENTRICITY
C            W=MASS OF ORBITED BODY (LB)
C            P=ORBITAL PERIOD (HRS)
C            TP=FRACTION OF ORBITAL PERIOD
C            THETA=ANGLE FROM PERIAPSIS (RADIANS)
C
      DIMENSION TP(100),THETA(100)
C
      DTH=PI(2.)/99.
      C=((1.-E*E)**1.5)/PI(2.)
       JRKS=1
       ST=0.0
       TH=0.0
      DO 200 J=1,100
       TP(J)=ST
       THETA(J)=TH
   10     F=1.+E*COS(TH)
       DST=C/(F*F)
       CALL RSFDX(ST,DST,DS,DTH,STS,JRKS)
       CALL RSCON(TH,DTH,JRKS)
      IF(JRKS.GT.1) GO TO 10
  200 CONTINUE
C
          F=1./TP(100.)
      DO 250 J=1,100
           TP(J)=TP(J)*F
  250 CONTINUE
C
      P=PI(2.)*(A**1.5)/(3600.*SQRT(1.068E-9*W))
                        RETURN
                         END
      SUBROUTINE SPHERE2(XO,PHI,RA,R,X,A,VN,PHIMIN,PHIMAX,N,NA)
C                    SUBROUTINE FOR GENERATION OF SHFAC
C                    PARAMETERS IN HELIOCENTRIC COORDINATES
C                    FOR A PORTION OF EARTH BETWEEN PHIMIN AND PHIMAX
C
      DIMENSION XO(3),X(3,520),A(520),VN(3,520),XC(3)
      DATA XC /3*0.0/
C
C                       CALCULATE SPHERICAL SURFACE TENSORS
C                       IN SPHERE'S COORDINATE SYSTEM
             CALL SPHERE(XC,R,X,A,VN,PHIMIN,PHIMAX,N,NA)
C
C                       TRANSFORM TO GENERAL COORDINATE SYSTEM
      DO 20 J=1,NA
C                        ROTATE ABOUT X AXIS
      CALL RTOP(X(2,J),X(3,J),RX,TD)
       RN=SQRT(VN(2,J)*VN(2,J)+VN(3,J)*VN(3,J))    $    TD=TD+PHI
      CALL PTOR(X(2,J),X(3,J),RX,TD)
      CALL PTOR(VN(2,J),VN(3,J),RN,TD)
C
C                        ROTATE ABOUT Z AXIS
      CALL RTOP(X(1,J),X(2,J),RX,TD)
       RN=SQRT(VN(1,J)*VN(1,J)+VN(2,J)*VN(2,J))    $    TD=TD+RA
      CALL PTOR(X(1,J),X(2,J),RX,TD)
      CALL PTOR(VN(1,J),VN(2,J),RN,TD)
C
C                        ROTATE FOR EARTH'S INCLINATION
      CALL RTOP(X(2,J),X(3,J),RX,TD)
      RN=SQRT(VN(2,J)*VN(2,J)+VN(3,J)*VN(3,J))    $    TD=TD-0.409274
      CALL PTOR(X(2,J),X(3,J),RX,TD)
      CALL PTOR(VN(2,J),VN(3,J),RN,TD)
C
C                        TRANSLATE
      X(1,J)=X(1,J)+XO(1)
      X(2,J)=X(2,J)+XO(2)
      X(3,J)=X(3,J)+XO(3)
C
   20 CONTINUE
            RETURN
             END
      REAL FUNCTION CPKAP(N)
C                  HEAT CAPACITY OF POLYIMIDE
C                       BTU/LB-R
      COMMON /FTAC/C(25,25),CAP(25),Q(25),T(25)
             CPKAP=4.6304E-4*T(N)
                      RETURN
                       END
      REAL FUNCTION CPGLAS(N)
C                  HEAT CAPACITY OF FUSED SILICA GLASS
C                        BTU/LB-HR
      COMMON /FTAC/C(25,25),CAP(25),Q(25),T(25)
      DIMENSION TR(5),CP(5)
C
      DATA TR/260.,860.,1260.,1960.,2460./
      DATA CP/0.103,0.23,0.261,0.285,0.285/, J/2/
C
           CPGLAS=YRDR(T(N),TR,CP,J,5)
                   RETURN
                    END
      REAL FUNCTION ETA(T,XLIFE)
C                       SOLAR CELL EFFICIENCY/TEMP FUNCTION
      DIMENSION TR(4),ETAR(4)
      DATA TR/366.,542.,600.,618./,   ETAR/0.16949,0.11637,0.09902,
     *                                     0.09362/,                J/2/
C
      ETA=(1.+0.036*(1.-XLIFE))*YRDR(T,TR,ETAR,J,4)
                       RETURN
                         END
      REAL FUNCTION FALBEDO(XA,AA,VNA,NA,XE,AE,VNE,FES,NE)
C               SUBROUTINE FOR SHAPE FACTOR  OF THE SUN
C               AS REFLECTED OFF THE EARTH
C
      DIMENSION XA(3,NA),AA(NA),VNA(3,NA), XE(3,NE),AE(NE),VNE(3,NE),
     *          XO(3),AO(1),VNO(3),FES(NE)
C
         F=0.0
      DO 10 J=1,NE
        IF(FES(J).EQ.0.0) GO TO 10
            XO(1)=XE(1,J)       $   VNO(1)=VNE(1,J)
            XO(2)=XE(2,J)       $   VNO(2)=VNE(2,J)
            XO(3)=XE(3,J)       $   VNO(3)=VNE(3,J)
                           AO(1)=AE(J)
      F=F+FES(J)*SHFAC(XA,AA,VNA,NA,XO,AO,VNO,1,XO,AO,VNO,1)
   10                      CONTINUE
               FALBEDO=F
                            RETURN
                             END
      REAL FUNCTION CRAD(E1,E2,A1,A2,F12)
                   C=0.0
                    IF(F12.GT.0.0)
     *C=1.714E-9*A1/(1./E1 -1.  +  1./F12  +  A1*(1./E2 -1.)/A2)
                CRAD=C
         RETURN
          END
      SUBROUTINE STDTM(DT,TIME,JRKS)
C                   SUBROUTINE FOR STABLE NUMERICAL INTEGRATION
C                   OF SOLAR ARRAY TEMPERATURE
C
      COMMON /FTAC/C(25,25),CAP(25),Q(25),T(25)
C
      C2=C(1,2)
      C4=C(1,4)
         T24=T(2)*T(2)*T(2)*T(2)
         T44=T(4)*T(4)*T(4)*T(4)
C
            TSS=SQRT(SQRT((Q(1)+C2*T24+C4*T44)/(C2+C4)))
            B=0.-(TSS*TSS+T(1)*T(1))*(TSS+T(1))*(C2+C4)/CAP(1)
C
      CALL TRSFDX(T(1),B,TSS,DTS,DT,TO,JRKS)
      CALL RSCON(TIME,DT,JRKS)
                    RETURN
                     END
      REAL FUNCTION SHFAC(X1,A1,VN1,N1,X2,A2,VN2,N2,XB,AB,VNB,NB)
C
C                   **SHAPE FACTOR FUNCTION**
C
C          REVISION OF 1/5/79
C              ... BLOCKAGE EVALUATED AFTER COSINE EVALUATION
C
C
C               X=COORDINATE ARRAY FOR EACH NODAL PLANE
C               A=AREA OF EACH NODAL PLANE
C               VN=COMPONENTS OF UNIT NORMAL VECTOR FOR EACH NODAL PLANE
C               N=NUMBER OF NODAL PLANES IN EACH BODY
C               1,2=BODIES FOR SHAPE FACTOR COMPUTATION
C               B=  ALL SURFACES BETWEEN BODY 1 AND BODY 2
C                 (IF NO INTERVENING SURFACES,LET XB=X1,ETC AND NB=1)
C               SHFAC=BLACK BODY SHAPE FACTOR REFERENCED TO BODY 1
C
      DIMENSION X1(3,N1),A1(N1),VN1(3,N1),X2(3,N2),A2(N2),VN2(3,N2),R(3)
      DIMENSION XB(3,NB),AB(NB),VNB(3,NB),B(3)
C
      A=0.0
      F1=0.0
C
C
      DO 10 J=1,N1
C
      A=A+A1(J)
C
      DO 20 K=1,N2
C               CALCULATE RADIUS VECTOR
      RMAG=0.0
      DO 30 LR=1,3
      R(LR)=X2(LR,K)-X1(LR,J)
   30 RMAG=RMAG+(R(LR)*R(LR))
      RMAG=SQRT(RMAG)
C
C
C               CALCULATE COSINES
      COS1=((VN1(1,J)*R(1))+(VN1(2,J)*R(2))+(VN1(3,J)*R(3)))/RMAG
      COS2=0.-(((VN2(1,K)*R(1))+(VN2(2,K)*R(2))+(VN2(3,K)*R(3)))/RMAG)
C
C               IF A COSINE IS .LE.0.0,THE SURFACES DON'T SEE
C               EACH OTHER
      IF(COS1.LE.0.0.OR.COS2.LE.0.0) GO TO 20
C
C                   BLOCKAGE EVALUATION
                          IF(NB.LE.1) GO TO 200
          DO 180 I=1,NB
C                           RADIUS VECTOR TO BLOCKING ELEMENT
           DO 150 LR=1,3
            B(LR)=XB(LR,I)-X1(LR,J)
  150      CONTINUE
C
           XDOTN=B(1)*VNB(1,I)+B(2)*VNB(2,I)+B(3)*VNB(3,I)
           RDOTN=R(1)*VNB(1,I)+R(2)*VNB(2,I)+R(3)*VNB(3,I)
                          IF(RDOTN.EQ.0.0) GO TO 180
           C=XDOTN/RDOTN
                          IF(C.LE.0.0.OR.C.GE.1.0) GO TO 180
C                   BLOCKAGE SURFACE VECTOR MAGNITUDE
                  S=0.0
                      DO 160 LR=1,3
                          S=S+(C*R(LR)-B(LR))**2
  160                 CONTINUE
                          IF(S.LE.AB(I)/2.) GO TO 20
  180 CONTINUE
  200 CONTINUE
C
      F1=F1+(COS1*COS2*A1(J)*A2(K)/(RMAG*RMAG))
   20 CONTINUE
   10 CONTINUE
C               SHAPE FACTOR
               SHFAC=F1/PI(A)
C
      RETURN
      END
      SUBROUTINE RTOP(X,Y,R,TH)
C                 RECTANGULAR TO POLAR CONVERSION
C
      R=SQRT(X*X+Y*Y)
               IF(R.GT.0.0) GO TO 10
                  TH=0.0
                  RETURN
   10 IF(X) 40,30,20
   20   IF(Y.LT.0.0) GO TO 25
             TH=ATAN(Y/X)
             RETURN
   25        TH=PI(2.)-ATAN(0.-Y/X)
             RETURN
   30        TH=PI(0.5)  $  IF(Y.LT.0.0)  TH=PI(1.5)
             RETURN
   40   IF(Y.LT.0.0) GO TO 45
             TH=PI(1.)+ATAN(Y/X)
             RETURN
   45        TH=PI(1.)-ATAN(0.-Y/X)
             RETURN
              END
      SUBROUTINE PTOR(X,Y,R,TH)
C                 POLAR TO RECTANGULAR CONVERSION
C
       X=R*COS(TH)
       Y=R*SIN(TH)
                     RETURN
                      END
      SUBROUTINE SPHERE(XO,R,X,A,VN,PHIMIN,PHIMAX,N,NA)
C
C                    SUBROUTINE FOR GENERATION OF SHFAC
C                    PARAMETERS FOR A SPHERE
C
      DIMENSION XO(3),X(3,520),A(520),VN(3,520)
C
C              XO=ORIGIN VECTOR FOR CENTER
C              R =RADIUS
C              X =SURFACE POSITION TENSOR
C              A =SURFACE AREA VECTOR
C              VN=SURFACE NORMAL TENSOR
C              PHIMIN,PHIMAX=LIMITS OF LATITUDE MEASURED FROM POLE (RADIANS)
C              N=NUMBER OF LATITUDE DIVISIONS
C              NA=ACTUAL NUMBER OF SURFACE ELEMENTS
C
    1 CONTINUE
      NA=0.0
      J=1
      PHI=PHIMIN
      DPHI=(PHIMAX-PHIMIN)/N
C
         THE=0.0
        DO 10 K=1,N
         SINP=SIN(PHI+DPHI/2.)
         COSP=COS(PHI+DPHI/2.)
          DA=DPHI/SINP
          XM=PI(2./DA)+0.5
          M=XM
          DA=PI(2./M)
         THE=THE+DA/2.
         AREA=R*R*DA*(COS(PHI)-COS(PHI+DPHI))
         NA=NA+M
                        IF(NA.LE.520) GO TO 15
C
                           I=(N-1)/2
                           N=2*I
                  WRITE(6,100) N
  100             FORMAT(1H0,20X,"SUBROUTINE SPHERE ..."/1H ,25X,
     *                   "CHOICE OF LAT. DIVISIONS PRODUCES EXCESS ELEME
     *NT ARRAY SIZE, N REDUCED TO",I4)
                  GO TO 1
C
   15 CONTINUE
          DO 20 L=1,M
           VN(1,J)=SINP*COS(THE+DA/2.)
           VN(2,J)=SINP*SIN(THE+DA/2.)
           VN(3,J)=COSP
           A(J)=AREA
           X(1,J)=XO(1)+R*VN(1,J)
           X(2,J)=XO(2)+R*VN(2,J)
           X(3,J)=XO(3)+R*VN(3,J)
              THE=THE+DA
              J=J+1
   20     CONTINUE
         PHI=PHI+DPHI
   10   CONTINUE
      RETURN
      END
      SUBROUTINE RSFDX(Y,DYDX,DY,DX,YS,KR)
C             ***RUNGE-KUTTA INTEGRATION SUBROUTINE***
      GO TO (10,12,12,13),KR
   10 YS=Y
      DY=DYDX*DX
      Y=YS+(DY/2.)
      GO TO 20
   12 DY=DY+(2.*DYDX*DX)
      IF(KR.LE.2) GO TO 14
      GO TO 15
   13 DY=(DY+(DYDX*DX))/6.
      Y=YS+DY
      GO TO 20
   14 Y=YS+(DYDX*DX/2.)
      GO TO 20
   15 Y=YS+(DYDX*DX)
   20 RETURN
      END
      SUBROUTINE TRSFDX(T,B,TSS,DTS,DELT,TO,KR)
C
C                *** TRANSFORM/RUNGE-KUTTA INTEGRATION ***
C                *** DEVELOPED BY  R.A.PROESCHEL- 8/77 ***
C                *** SUCESSIVE APPROXIMATION VERSION OF 4/78 ***
C
                      GO TO (1,2,3,4),KR
C
    1 DTS=TDELT(1.,0.,B,DELT)
      TO=T
      T=TDELT(TSS,TO,B,DELT/2.)+TO
        RETURN
    2 DTS=TDELT(1.,0.,B,DELT)*2.+DTS
      T=TDELT(TSS,TO,B,DELT/2.)+TO
        RETURN
    3 DT =TDELT(1.,0.,B,DELT)
      DTS=DTS+2.*DT
      T=TO+(TSS-TO)*DT
        RETURN
    4 T=TO+(TSS-TO)*(DTS+TDELT(1.,0.,B,DELT))/6.
        RETURN
         END
      REAL FUNCTION TDELT(TSS,TO,B,DELT)
       PHI=B*DELT
        IF(PHI.LT.-25.) GO TO 10
         TD=(TSS-TO)*(1.-EXP(PHI))
        GO TO 20
   10    TD=TSS-TO
   20 TDELT=TD
      RETURN
      END
      SUBROUTINE RSCON(XC,DXC,KC)
C             ***CONTROL DECK FOR RUNGEKUTTA INTEGRATION***
      GO TO (10,12,13,14),KC
   10 XC=XC+(DXC/2.)
      KC=2
      GO TO 20
   12 KC=3
      GO TO 20
   13 KC=4
      XC=XC+(DXC/2.)
      GO TO 20
   14 KC=1
   20 RETURN
      END
      REAL FUNCTION YRDR(X1,XT,YT,J,N)
C        ***REVISED HIGH SPEED VERSION***
      DIMENSION XT(N),YT(N)
   50 IF(X1.GT.XT(J)) GO TO 30
      IF(X1.GT.XT(J-1).OR.J.LE.2) GO TO 60
      J=J-1
      GO TO 50
   30 J=J+1
      IF(J.LE.N) GO TO 50
      J=J-1
   60 YRDR=YT(J-1)+(((YT(J)-YT(J-1))/(XT(J)-XT(J-1)))*(X1-XT(J-1)))
      RETURN
      END
      REAL FUNCTION PI(X)
C                 ** FUNCTION SIMULATING ROCKWELL SCCLIB **
C                 **              ROUTINE "PI"           **
      PI=3.1415927*X
       RETURN
        END
